{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Risk Premia Estimation using GMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by importing the modules and functions needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import hstack, ones, array, mat, tile, reshape, squeeze, eye, asmatrix\n",
    "from numpy.linalg import inv\n",
    "from pandas import read_csv, Series \n",
    "from scipy.linalg import kron\n",
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next a callable function is used to produce iteration-by-iteration output when using the non-linear optimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iteration = 0\n",
    "last_value = 0\n",
    "function_count = 0\n",
    "\n",
    "def iter_print(params):\n",
    "    global iteration, last_value, function_count\n",
    "    iteration += 1\n",
    "    print(f'Func value: {last_value:6.6g}, Iteration: {iteration}, Function Count: {function_count}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GMM objective, which is minimized, is defined next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gmm_objective(params, p_rets, f_rets, Winv, out=False):\n",
    "    global last_value, function_count\n",
    "    t,n = p_rets.shape\n",
    "    t,k = f_rets.shape\n",
    "    beta = squeeze(array(params[:(n*k)]))\n",
    "    lam = squeeze(array(params[(n*k):]))\n",
    "    beta = reshape(beta,(n,k))\n",
    "    lam = reshape(lam,(k,1))\n",
    "    betalam = beta @ lam\n",
    "    expected_ret = f_rets @ beta.T\n",
    "    e = p_rets - expected_ret\n",
    "    instr = tile(f_rets,n)\n",
    "    moments1  = kron(e,ones((1,k)))\n",
    "    moments1 = moments1 * instr\n",
    "    moments2 = p_rets - betalam.T\n",
    "    moments = hstack((moments1,moments2))\n",
    "\n",
    "    avg_moment = moments.mean(axis=0)\n",
    "    \n",
    "    J = t * mat(avg_moment) * mat(Winv) * mat(avg_moment).T\n",
    "    J = J[0,0]\n",
    "    last_value = J\n",
    "    function_count += 1\n",
    "    if not out:\n",
    "        return J\n",
    "    else:\n",
    "        return J, moments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `G` matrix, which is the derivative of the GMM moments with respect to the parameters, is defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gmm_G(params, p_rets, f_rets):\n",
    "    t,n = p_rets.shape\n",
    "    t,k = f_rets.shape\n",
    "    beta = squeeze(array(params[:(n*k)]))\n",
    "    lam = squeeze(array(params[(n*k):]))\n",
    "    beta = reshape(beta,(n,k))\n",
    "    lam = reshape(lam,(k,1))\n",
    "    G = np.zeros((n*k+k,n*k+n))\n",
    "    ffp = (f_rets.T @ f_rets) / t\n",
    "    G[:(n*k),:(n*k)]=kron(eye(n),ffp)\n",
    "    G[:(n*k),(n*k):] = kron(eye(n),-lam)\n",
    "    G[(n*k):,(n*k):] = -beta.T\n",
    "    \n",
    "    return G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the data is imported and a subset of the test portfolios is selected to make the estimation faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = read_csv('FamaFrench.csv')\n",
    "\n",
    "# Split using both named colums and ix for larger blocks\n",
    "dates = data['date'].values\n",
    "factors = data[['VWMe','SMB','HML']].values\n",
    "riskfree = data['RF'].values\n",
    "portfolios = data.iloc[:,5:].values\n",
    "\n",
    "t,n = portfolios.shape\n",
    "portfolios = portfolios[:,np.arange(0,n,2)]\n",
    "t,n = portfolios.shape\n",
    "excess_ret = portfolios - np.reshape(riskfree,(t,1))\n",
    "k = np.size(factors,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting values for the factor loadings and rick premia are estimated using OLS and simple means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betas = []\n",
    "for i in range(n):\n",
    "    res = sm.OLS(excess_ret[:,i],sm.add_constant(factors)).fit()\n",
    "    betas.append(res.params[1:])\n",
    "\n",
    "avg_return = excess_ret.mean(axis=0)\n",
    "avg_return.shape = n,1\n",
    "betas = array(betas)\n",
    "res = sm.OLS(avg_return, betas).fit()\n",
    "risk_premia = res.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The starting values are computed the first step estimates are found using the non-linear optimizer.  The initial weighting matrix is just the identify matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risk_premia.shape = 3\n",
    "starting_vals = np.concatenate((betas.flatten(),risk_premia))\n",
    "\n",
    "Winv = np.eye(n*(k+1))\n",
    "args = (excess_ret, factors, Winv)\n",
    "iteration = 0\n",
    "function_count = 0\n",
    "opt = minimize(gmm_objective, starting_vals, args=args, callback=iter_print)\n",
    "step1opt = opt.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we look at the risk premia estimates from the first step (inefficient) estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "premia = step1opt[-3:]\n",
    "premia = Series(premia,index=['VWMe', 'SMB', 'HML'])\n",
    "print('Annualized Risk Premia (First step)')\n",
    "print(12 * premia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the first step estimates are used to estimate the moment conditions which are in-turn used to estimate the optimal weighting matrix for the moment conditions.  This is then used as an input for the 2nd-step estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = gmm_objective(step1opt, excess_ret, factors, Winv, out=True)\n",
    "S = np.cov(out[1].T)\n",
    "Winv2 = inv(S)\n",
    "args = (excess_ret, factors, Winv2)\n",
    "\n",
    "iteration = 0\n",
    "function_count = 0\n",
    "opt = minimize(gmm_objective, step1opt, args=args, callback=iter_print)\n",
    "step2opt = opt.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally the VCV of the parameter estimates is computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = gmm_objective(step2opt, excess_ret, factors, Winv2, out=True)\n",
    "G = gmm_G(step2opt, excess_ret, factors)\n",
    "S = np.cov(out[1].T)\n",
    "vcv = inv(G @ inv(S) @ G.T)/t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The annualized risk premia and their associated t-stats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "premia = step2opt[-3:]\n",
    "premia = Series(premia,index=['VWMe', 'SMB', 'HML'])\n",
    "premia_vcv = vcv[-3:,-3:]\n",
    "print('Annualized Risk Premia')\n",
    "print(12 * premia)\n",
    "\n",
    "premia_stderr = np.diag(premia_vcv)\n",
    "premia_stderr = Series(premia_stderr,index=['VWMe', 'SMB', 'HML'])\n",
    "print('T-stats')\n",
    "print(premia / premia_stderr)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
